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SUMMARY 


A description of an improved version of the NASA/Lockheed multielement airfoil 
analysis computer program is presented. The improvements include several major 
modifications of the aerodynamic model as well as substantial changes of the 
computer code. The modifications of the aerodynamic model comprise the repre- 
sentation of the boundary layer and wake displacement effects with an equiva- 
lent source distribution, the prediction of wake parameters with Green's lag- 
entrainment method, the calculation of turbulent boundary layer separation with 
the method of Nash and Hicks, the estimation of the onset of confluent boundary 
layer separation with a modified form of Goradia's method, and the prediction 
of profile drag with the formula of Squire and Young. The paper further de- 
scribes the modifications of the computer program for which the structured ap- 
proach to computer software development was employed. Important aspects of the 
structured program development such as the functional decomposition of the 
aerodynamic theory and its numerical implementation, the analysis of the data 
flow within the code, and the application of a pseudo code are discussed. 

Computed results of the new program version are compared with recent experi- 
mental airfoil data. The comparisons include global airfoil parameters such as 
lift, pitching moment and drag coefficients, and distributions of surface pres- 
sures and boundary layer velocity profiles. 


INTRODUCTION 


In the past, high lift design and technology rested in the hands of a few ex- 
perienced aerodynami cists. Design methodology and criteria were heavily in- 
fluenced by the analytical inviscid flow methods and the experimental data a- 
vailable. With the advent of high-speed computers and the appearance of im- 
proved models for turbulent flows, many complex problems, including high-lift 
design and analysis, were attacked theoretically. 


The work reported in this paper was supported partly by NASA-Langley con 
tract NASI -14522 and partly by the Independent Research and Development 
Program of The Boeing Company. 
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I 


One such approach to high-lift, multielement airfoil analysis was developed at 
Lockheed-Georgia under the sponsorship of the NASA-Langley Research Center 
(ref. 1). This program was among the first attempts at analyzing the complex 
viscous flow about slotted airfoils and has received worldwide distribution and 
usage. A unique feature of this multielement airfoil program is the model of 
the confluent boundary layer flow (ref. 2). 

Over the years, the original version of the program was modified extensively to 
improve its predictions for different types of high-lift airfoils. Many im- 
provements, mainly in the area of the potential flow calculation, were made by 
researchers at the Langley Research Center (ref. 3). For this reason, the code 
is generally referred to as the NASA/Lockheed multielement airfoil program. A 
version for single element airfoils was recently extracted from the multiele- 
ment airfoil code by researchers at North Carolina State University (ref. 4). 

Widespread and steady usage of the computer program clarified its strengths and 
weaknesses. Both favorable and unfavorable aspects have been brought to the 
surface by continued attempts at using the program as an engineering tool. The 
more serious shortcomings were the lack of agreement between the documentation 
and the available version of the code and the high failure rate in applying the 
method for various configurations. However, the program was found to contain 
sufficient positive features to justify its choice as a starting point for ad- 
ditional theoretical work in the high-lift area. 

This paper briefly describes the aerodynamic theory and the corresponding com- 
puter program of a new version of the multielement airfoil program; a detailed 
description can be found in references 5 and 6. Symbols are defined in an 
appendix. 

MULTIELEMENT AIRFOILS 


The flow around high-lift airfoils is characterized by many different inviscid 
and viscous flow regions. Their complex physics is illustrated by figure 1. 

In particular, the existence of confluent boundary layers and the regions of 
separated flow distinguish the high-lift airfoil problem from the aerodynamic 
problem of airfoils at cruise conditions. The various flow regions, including 
the outer potential flow, the ordinary laminar and turbulent boundary layers, 
viscous wakes, and the confluent boundary layer, are analyzed by the code. 
Furthermore, the prediction of transition from laminar to turbulent boundary 
layer flow and the prediction of the onset of boundary layer separation are a 
necessary part of the code. Cove separation and large scale separation phe- 
nomena, however, are not modeled. 


PROGRAM MODIFICATIONS 


The new program version differs from the baseline version (ref. 3) in the fol- 
lowing areas: 
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1) The method used to represent the effect of the viscous flow on the outer 
potential flow, termed equivalent airfoil representation in the baseline 
version of the program, has been modified. It has been •replaced by the 
surface transpiration method which uses a distribution of sources along 
airfoil surfaces and wake centerlines to model boundary layer and wake 
displacement effects. 

2) The flow model of the potential core region has been changed. The new 
method performs independent boundary layer and wake calculations. These 
calculations utilize the ordinary laminar and turbulent boundary layer 
routines of the baseline version of the code, and in addition, the lag-en- 
trainment method of reference 7 for wake flows. The revised flow model 

of the core region calculates the location of the wake centerlines. 

3) An attempt is made to predict the onset of separation of the confluent 
boundary layer by a modified version of Goradia's confluent boundary layer 
method. In this method, the power law velocity profile of the wall layer 
is replaced by Coles' two-parameter velocity profile (ref. 8). 

4) The drag prediction method of Squire and Young (ref. 9) has been incorpo- 
rated into he program, replacing the previous pressure and skin friction 
integration scheme. 

5) The original method used for the prediction of separation for ordinary 
turbulent boundary layer flow has been replaced b> the Boeing version of 
the method of Nash and Hicks (ref. 10). 

6) The modifications of the aerodynamic theory required a major overhaul of 
the computer code. Most parts of the code have been rewritten using a 
systematic approach to computer software design. This work was guided by 
a functional decomposition of the many aspects of the aerodynamic model 
and its numerical implementation. In addition, a detailed study was made 
of the data flow within the program, and the logic of the code was out- 
lined prior to the actual program development using a pseudo code. The 
most important aspects of this work are briefly reviewed in this paper. 


AERODYNAMIC FLOW MODELS 


The aerodynamic theory of the new version of the computer program is outlined 
below with emphasis on modifications. The aerodynamic analysis and its numer- 
ical implementation assume two-dimensional, subsonic flow in which all boundary 
layers are attached to the airfoil surface. 


Potential Flow 


Inviscid, irrotational flow is calculated using the stream function approach of 
Oeller (ref. 11). Laplace's equation is solved subject to the boundary condi- 
tion of a constant value of the stream function on airfoil surfaces. The meth- 
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od is of the panel type, see figure 2, with a constant strength vortex distri- 
bution on each panel of the airfoil surface. A formulation of the Kutta condi- 
tion is imposed which requires the tangential velocities at the upper and lower 
surface trailing edge points to be equ?l. Compressibility effects are taken 
into account by employing the Karman-Tsien rule. 


Viscous Flow Representation 


Oeller's stream function method has been modified in order to account for dis- 
placement effects of boundary layers and wakes within the potential flow solu- 
tion. An equivalent distribution of sources simulating the viscous flow dis- 
placement thickness is placed on the surface and wake centerline of an airfoil 
component. This is the surface transpiration method which, within the frame- 
work of thin boundary layer theory, is completely equivalent to the method of 
geometrically adding the displacement thickness to the basic airfoil geometry. 
Application of this technique is computationally efficient, since most of the 
aerodynamic influence coefficients do not c. ige during the solution procedure. 

The strength a of the equivalent source distribution is obtained from 



where 6* denotes the displacement thickness of either boundary layer or wake, 
and the symbol U stands for the inviscid flow velocity on airfoil surface and 
wake centerline. The variable s represents arc length. The computed source 
distribution is discretized using panels with constant source strength on the 
surface and wake centerline of each airfoil component. 

It should be emphasized that the employed flow model does not account for wake 
curvature e-ffects. Consequently, constant strength vortex panels are only used 
on the surface of an airfoil and not on its wake centerline, see figure 2. 


Wake Centerline 


The capability of computing the position of wnke centerlines has been added to 
the program as part of the revision of the flow model in the core region, see 
figure 1. A wake centerline is part of the stagnation streamline. 

Since the potential flow problem is solved on the basis of a stream function 
approach, which in addition to the surface velocity provides the value of the 
stream function for each stagnation streamline, it is convenient to also use 
the stream function formulation to trace wake centerlines. This is done in an 
iterative procedure beginning with an assumed initial position of a wake cen- 
terline. During each step of this iteration the locations of all panels of the 
wake centerline are updated simultaneously by solving a linearized form of the 
following stream function equation: 
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'Pia = P (x) 

Here, ’I'm denotes the known value of the stream function at a stagnation stream- 
line. The variable x represents the array of unknown panel corner point co- 
ordinates of the wake centerline. 


Laminar Boundary Layer 


Laminar boundary layer characteristics are calculated with the compressible 
method of Cohen and Reshotko (ref. 12) who reduced the problem to simple quad- 
rature. Application of a compressible method seems to be necessary for slot- 
ted high-lift airfoils, since laminar boundary layers often exist in the slot 
between neighboring airfoil compor.ants, where even at low free stream Mach num- 
bers the flow is highly compressible with velocities approaching and frequently 
exceeding sonic conditions. 

Laminar separation is predicted with the criterion of Goradia and Lyman (ref. 
13) which is an empirical correlation of the local values of the Mach number 
gradient with the momentum thickness Reynolds number. Either laminar stall or 
the occurrence of laminar short bubble separation is predicted. In the case 
of laminar short bubbles, subsequent turbulent reattachment of the separated 
laminar boundary layer is assumed, but neither the length of the separation 
bubble nor the details of the flow within the bubble are modeled. 

The computer program provides two options for transition from laminar to tur- 
bulent boundary layer flow. The user can eithe* specify fixed transition 
points, such as the location of trip strips, or can compute free transition. 

In the latter case, a standard two-step approach is employed. 


Turbulent Boundary Layer 


Two different integral methods determine the characteristics of ordinary tur- 
bulent boundary layers. The method of Truckenbrodt (ref. 14) is used during 
the iterative solution procedure. It is an incompressible approach based on 
the momentum and energy integral equations. Goradia (ref. 1) introduced the 
idea of constraining the shape factor H, defined as the ratio of energy dis- 
sipation thickness to momentum thickness, in order to avoid premature separa- 
tion of the turbulent boundary layer during the first cycles of the iteration. 
This approach avoids failures of the boundary layer integration at separation 
and can be viewed as an artificial way of modeling separated flows. The inte- 
gral method of Nash and Hicks (ref. 10), which accounts for the history of tur- 
bulent s..aar stresses, is applied at the end of the iterative solution proce- 
dure for the purpose of computing boundary layer separation. Displacement 
thickness and skin friction obtained from the method of Nash and Hicks are not 
utilized. 
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Wake Flow 


The properties of turbulent wakes are analyzed with the lag-entrainment method 
of reference 7. The method is formulated in terms of the momentum integral 
equation, the entrainment equation, and an empirical equation for the *.tream- 
wise rate of change of the entrainment coefficient. The entrainment equation 
is derived from the definition of the entrainment coefficient, which represents 
the change of mass flow within the wake layer. An incompressible version of 
Green's treatment of wake flow is used, neglecting the effects of curvature on 
the mean flew and the turbulence structure of the wake. 


Confluent Boundary Layer 


The program computes confluent boundary layers with the model of Goradia (ref. 
2). In this model, the confluent boundary layer downstream of the core region 
is divided into two regions. In the first region, turbulent mixing of wake and 
boundary layer is incomplete. The mean velocity profile clearly shows the re- 
mainder of the wake profile, see figure 3. In the second region, the effect 
of the wake is not visible in the mean velocity profile which is similar to 
that of a will jet. Downstream of the second region, the confluent boundary 
layer degenerates into an ordinary turbulent boundary layer. 

Goradia formulated an integral method by subdividing the confluent boundary lay- 
er into several layers and assuming the validity of boundary layer equations 
and self-similarity of the mean velocity profile in each of these layers. The 
method is incompressible, neglects curvature effects, and relies to a large ex- 
tent on empirical information about shear stress, the rate of growth of vari- 
ous layers, and velocity profiles. Furthermore, the ntodel ignores multiple 
wakes and multiple potential cores that might exist near the trailing edgp of a 
high-lift airfoil consisting of more than two components. The shape factor H 
of the layer adjacent to the airfoil surface, termed the wall layer, is con- 
strained in order to avoid program failures in regions of separated flow. For 
this reason, Goradia's confluent boundary layer method is applied during the i- 
teration procedure, when unrealistic potential flow pressure distributions can 
cause premature boundary layer separation. 

The described flow model has been modified to predict separation of confluent 
boundary layers. The power law velocity profile of the wall layer has been ; e- 
placed by Coles' tv<o-parameter profile (ref. 8), which is known to provide a 
realistic representation of ordinary turbulent boundary layers near separation. 
The shape factor of the wall layer is not constrained in this modification, but 
most other 'eatures of Goradia's confluent boundary layer model including its 
empirical content are retained. 


iterative Solution Procedure 


The solution is iterative, since most of the individual flow problems and their 
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coupling are nonlinear. The computer program uses a conventional cyclic itera- 
tion procedure in which each cycle consists of the following steps: 

Step 1 : A potential flow solution for the multielement airfoil is calcu- 

lated. During the first cycle of the iteration, the calculation 
is performed without any representation of viscous flow displace- 
ment effects. 

Step 2: The positions of the wake centerlines are computed. 

Step 3: Solutions of all viscous flow problems including laminar and tur- 

bulent boundary layers, confluent boundary layers, and viscous 
wakes are calculated with the potential flow velocities and wake 
centerline locations obtained in the previous steps as input data. 

At the end of this computational step, the displacement thicknesses 
of all boundary layers and wakes are available. 

Step 4: A source distribution representing the displacement effect of all 

viscous layers is computed. 

The computer program does not rely on a convergence criterion. Instead, five 
iteration cycles are always executed and the user of the program must judge the 
quality of the solution. 

In order to assist the iteration scheme in arriving at a converged solution, 
the source strength a computed in Step 4 of each cycle is modified by adding 
2/3 of a computed in this iteration cycle to 1/3 of a computed in the pre- 
vious iteration cycle. 


Forces and Moments 


At the end of each iteration cycle, airfoil lift and pitching moment coeffi- 
cients are computed by an integration of surface pressure and skin friction. 
Profile drag is obtained by applying the formula of Squire and Young (ref. 9) 
The total profile drag of a multielement airfoil is assumed to be the sum of 
the contributions of its components. The drag of each airfoil component in 
turn is calculated from the values of boundary layer momentum thickness and 
surface velocity at the component trailing edge. 


COMPUTER CODE 


The n®w version of the code is written in the CDC FORTRAN extended 4 (FTN4) 
language and will run under the CDC Network Operating System (NOS). The CDC 
overlay system is used to assure that the cide will execute in a field length 
less than 100 K octal . 

The programming methodologies used to design and develop the new version of the 
computer code include a functional d., :omposition of the aerodynamic theory, ^ 
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data flow analysis, and a control flow analysis. 

Each of these related design tasks was performed several times in an iterative 
manner to produce a final design for the new version of the computer code be- 
fore changes or improvements to the baseline code were made. The final design 
resulted in major changes of the baseline version of the program in the follow- 
ing sections: upper level control routines, geometry preprocessing routines, 

and the potential flow solution routines. The final design was also used to 
integrate the new aerodynamic models into the baseline code. Table 1 lists all 
subroutines in the baseline version of the code and indicates the type of 
changes made to incorporate them in the new version. 

The functional decomposition of the code was based on the engineering specifi- 
cation of the aerodynamic models and the numerical techniques necessary for 
their solution. Figure 4 shows the upper level decomposition chart where the 
major functions are defined in engineering terms. The complex physics of the 
flow about multielement airfoils is reflected in these charts. 

The data flow analysis of the code was done for each module identified in the 
functional decomposition by specifying the input and output data for the module 
as well as its own decomposition. Control of the data flow within a module is 
maintained by requiring that the input to any of its submodules must be either 
an input to the module or the output of another of its submodules. 

The control flow analysis of the new code was done with the aid of a pseudo 
code, which is a small set of simple logic and loop statements which suffice 
to describe the control within a module of the functional decomposition. Al- 
though the submodules of a module can be used in any sequence and any number of 
times to complete the function of the module, it is an aim of the design pro- 
cess to keep the control within a module as simple as possible. All new sub- 
routines in the code include as comment cards the pseudo code for the module 
which they implement. 


TEST-THEORY COMPARISONS 


In the following comparison of theoretical and experimental airfoil data, three 
versions of the NASA/Lockheed multielement airfoil program are referred to: 

Version A: This is the baseline version of the computer program with minor 

modifications. The baseline version was available from the NASA 
in June 1976. 


Version B: This version, described in reference 15, differs from version A 

in two areas. Profile drag is predicted by the Squire and Young 
formula. Separation of ordinary turbulent boundary layers is 
calculated using the method of Nash and Hicks. 

Version C: This is the version of the program described in this paper. 


170 


A large number of airfoil configurations were analyzed using the program ver- 
sions listed above. Only a few results of this program evaluation, concerning 
the GA(W)-1 airfoil and a Boeing multielement high-lift airfoil (figure 5), are 
discussed in this paper. A detailed report of this evaluation is contained in 
reference 6. 


6A(W)-1 Airfoil 


This airfoil was chosen to test the program capability of predicting perfor- 
mance characteristics of single airfoils. Figure 6 contains the theoretical 
lift, pitching moment, and drag curves and their comparison with experimental 
data of McGhee and Beasley (ref. 16). ?oth Version A and the new program Ver- 
sion C predict identical lift and moment curves that in turn agree well with 
measured GA(W)-1 data up to the onset of trailing edge stall at about 8 degrees 
angle of attack. Trailing edge stall is not modeled by any of the program ver- 
sions. 

Considerable differences between all drag polars are shown in figure 5. Ver- 
sion A, utilizing an integration of surface pressure and skin friction in the 
prediction of profile drag, gives the highest drag coefficients. Version C, 
applying the Squire and Young formula, offers drag values that are lower than 
the corresponding experimental drag coefficients. The lack of agreement of the 
three drag polars emphasizes the fact that even for single airfoils at low 
speed the problem of obtaining accurate drag computations is not yet solved. 

Surface pressures of the GA(W)-1 at 8 degrees angle of attack, computed by pro- 
gram Version C and plotted in figure 7, agree well with their experimental 
counterparts. In this figure, the symbols S and LS refer to theoretical points 
of turbulent separation and laminar short bubbles, respectively. The symbol FT 
indicates the experimental trip strip location which is specified as a fixed 
transition point in the computer simulation. A laminar short bubble with sub- 
sequent turbulent reattachment of the boundary layer is indicated near the up- 
per surface leading edge, and turbulent boundary layer separation is predicted 
theoretically near the upper surface trailing edge. The latter prediction is 
confirmed by the* experimental pressure distribution which shows a constant 
pressure downstream of the theoretical point of separation. 


Boeing High-Lift Airfoil 


The Boeing four-element high-lift airfoil was used as the main test case for 
multiple airfoils. It consists of a wing section with a leading edge flap and 
a double -slotted trailing edge flap. Global airfoil parameters and detailed 
distributions of surface pressures and boundary layer data are available for 
comparisons. The data were obtained in the Boeing Research Wind Tunnel (BRWT) 
on a model with a 2 foot unextended wing chord and a 5 foot span. Careful 
blowing of the wall boundary layers was applied in order to closely approximate 
a two-dimensional flow pattern across the wnole span of the airfoil. 
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The lift curve, pitching ir.O(nent, and drag polar of this airfoil at a Reynolds 
number of two million, based on the wing reference chord, are given in figure 
8. The experimental lift coefficients are balance data that are within 1.5% 
of the lift obtained by pressure integration. The profile drag of the airfoil 
is the result of wake rake measurements taken at a fixed spanwise position rel- 
atively free from interference effects of flap supporting brackets and pres- 
sures taps. The maximum spanwise variation of the measured airfoil drag is in- 
dicated in figure 8. 

All attempts failed when using program Version A to obtain a converged solution 
for this airfoil. Program Version B arrived at converged solutions between 8 
and 20 degrees angle of attack, but underpredicted the lift by a considerable 
amount. The theoretical predictions of program Version C match the experi- 
mental lift coefficients well at angles of attack below the onset of trailing 
edge stall, which is theoretically predicted by the program to take place at 
about 16 degrees. 

The theoretical values of the profile drag of Version C are relatively close to 
the measured profile drag. In judging the quality of this drag comparison, the 
reader should recall the problems of two-dimensional high-lift testing and the 
uncertainties In applying the Squire and Young formula to multielement air- 
foils. 

Figure 9 contains comparisons of theoretical and experimental surface pressures 
of wing and main flap at 8.4 degrees angle of attack. This figure confirms the 
earlier finding that Version C indeed provides the best theoretical results. 
Differences between the theory of Version C and experiment are noted in cove 
regions and on the upper surface of the wing near the leading edge. The latter 
problem is due to the failure of the program to accurately simulate the flow on 
the lower surface of the leading edge device. 

Figure 10 shows boundary layer velocity profiles at several chordwise stations 
on the upper wing surface. The experimental velocity profiles reveal that very 
little confluence of the wake behind the leading edge device and the wing 
bounc'ary layer has taken place and that an initially existing weak confluent 
boundary layer above the wing has degenerated early into an ordinary turbulent 
boundary layer. This feature of the flow field is very well simulated by Ver- 
sion (’, but not by Versions A and B. 


CONCLUSIONS 


The following conclusions about the reliability and quality of the predictions 
of the new program are drawn: 

1) The reliability of the program executions has been greatly improved. All 
test cases have produced converged solutions within a few iteration cycles. 
This improvement is a consequence of the applicatirn of the structured 
approach to computer programming where much attention was paid to the func- 
tional decomposition of the aerodynamic model, its numerical implementa- 
tion, and the data flow within the code. 
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The accuracy of the program predictions has been improved. This is due to 
several major modifications of the aerodynamic model - above a11, due to 
the different representation of the viscous flow displacement effects and 
the improved model of the potential core region. 

The computed results are consistent with the basic assumptions of the aero- 
dynamic model. Best results are obtained in cases where most of the flow 
is attached to the airfoil's surface, but the quality of the predictions 
gradually deteriorates with increasing trailing edge stall and cove separa- 
tion. 

The usefulness of the confluent boundary layer method of Goradia and its 
modification developed for the purpose of predicting confluent boundary 
layer separation have not yet been tested. Configurations were chosen for 
most of the program evaluation with little confluence of wakes and bound- 
ary layers. 

The performance of the program needs to be tested for configurations at off 
optimum shape design. 

The evaluation of the computer program was hampered by the shortage of re- 
liable experimental high-lift data. Additional wind tunnel testing of some 
of the more important high-lift airfoil configurations will Increase con- 
fidence in their predicted performance. 


APPENDIX 


SYMBOLS 


C airfoil reference chord 

Cj drag coefficient 

Cl lift coefficient 

Cju pitching moment coefficient about the quarter chord point 

Cp surface pressure coefficient 

H ratio of energy dissipation thickness and momentum thickness 

M* free stream Mach number 

R|^ Reynolds number formed by free stream velocity and airfoil reference 

chord 

s arc 1 ength 

U inviscid surface velocity or wake centerline velocity 

U* free st'"?am velocity 

u boundary layer velocity parallel to airfoil surface 

X X - coordinate of global axis system 

? array of panel corner points 

Y coordinate normal to airfoil surface 

a angle of attack 

6* displacement thickness 

a source strength 

stream function 

'I'm stream function value at a stagnation streamline 

Abbreviations 

FT fixed transition 

LS laminar short bubble 

S separation 
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Out«r Potential Flow 



Layer 

Figure 1.- Flow regions of multielement airfoils. 



Figure 2,- Potential-flow singularities. 






































